Reconciling Semiclassical and Bohmian Mechanics: 
II. Scattering states for discontinuous potentials 
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In a previous paper [J. Chem. Phys. 121 4501 (2004)] a unique bipolar decomposition, ^ = + 
5*2 was presented for stationary bound states 'I' of the one-dimensional Schrodinger equation, such 
that the components ^'i and ^2 approach their semiclassical WKB analogs in the large action limit. 
Moreover, by applying the Madelung-Bohm ansatz to the components rather than to ^ itself, the 
resultant bipolar Bohmian mechanical formulation satisfies the correspondence principle. As a result, 
the bipolar quantum trajectories are classical-like and well-behaved, even when ^I" has many nodes, or 
is wildly oscillatory. In this paper, the previous decomposition scheme is modified in order to achieve 
the same desirable properties for stationary scattering states. Discontinuous potential systems are 
considered (hard wall, step, square barrier/well), for which the bipolar quantum potential is found to 
be zero everywhere, except at the discontinuities. This approach leads to an exact numerical method 
for computing stationary scattering states of any desired boundary conditions, and reflection and 
transmission probabilities. The continuous potential case will be considered in a future publication. 



I. INTRODUCTION 

Much attention has been directed by theoreti- 
cal/computational chemists towards developing reliable 
and accurate means for solving dynamical quantum 
mechanics problems — i.e., for obtaining solutions to 
the time-dependent Schrodinger equation — for molecu- 
lar systems. Insofar as "exact" quantum methods are 
concerned, two traditional approaches have been used: 
(1) representation of the system Hamiltonian in a fi- 
nite, direct-product basis set; (2) discretization of the 
wavefunction onto a rectilinear grid of lattice points over 
the relevant region of configuration space. Both ap- 
proaches, however, suffer from the drawback that the 
computational effort scales exponentially with system 
dimensionality.^'^ Recently, a number of promising new 
methods have emerged with the potential to alleviate 
the exponential scaling problem once and for all. These 
include various basis set optimization methodsf^i^i^ 
and build-and-prune methods^ such as those based on 
wavelet techniques.— ^^ii^ 

On the other hand, a completely different approach 
to the exponential scaling problem is to use basis sets or 
grid points, that themselves evolve over time. The idea is 
that at any given point in time, one need sample a much 
smaller Hilbert subspace, or configuration space region, 
than would be required at all times — thus substantially 
reducing the size of the calculation. For basis set calcula- 
tions, much progress along these lines has been achieved 
by the multi-configurational time-dependent Hartree 
(MCTDH) method, developed by Meyer, Manthe and co- 
workers Jii^ More recently, time-evolving grid, or "quan- 
tum trajectory" methods^iiiii^ii^iilii^ (QTMs) have 
also been developed, and for certain types of systems, 
successfully applied at quite high dimensionalities .ililS 

QTMs are based on the hydrodynamical picture of 
quantum mechanics, developed over half a century ago 
by Bohm^^'^° and Takabayasi,^ who built on the earlier 



work of Madelung^^ and van Vlecki^ QTMs are inher- 
ently appealing for a number of reasons. First, they offer 
an intuitive, classical-like understanding of the underly- 
ing dynamics, which is difficult-to-impossible to extract 
from more traditional fixed grid/basis methods. In ef- 
fect, quantum trajectories are like ordinary classical tra- 
jectories, except that they evolve under a modified po- 
tential V + Q, where Q is the wavefunction-dependent 
"quantum potential" correction. Second, QTMs hold 
the promise of delivering exact quantum mechanical re- 
sults without exponential scaling in computational effort. 
Third, they provide a pedagogical understanding of en- 
tirely quantum mechanical effects such as tunnelin g^^i^^ 
and interference They have already been used to 
solve a variety of different types of problems, including 
barrier transmission^ non-adiabatic dynamics;^ and 
mode relaxationiSI Several intriguing phase space gener- 
alizations have also emerged j^^'^^i^^'^" of particular rele- 
vance for dissipative systems ^ii^i^^i^ 

Despite this success, QTMs suffer from a significant 
numerical drawback, which, to date, precludes a com- 
pletely robust application of these methods. Namely: 
QTMs are numerically unstable in the vicinity of am- 
plitude nodes. This "node problem" manifests in several 
different waysi ^^'^^ (1) infinite forces, giving rise to kinky, 
erratic trajectories; (2) compression/inflation of trajecto- 
ries near wavefunction local extrema/nodes, leading to; 
(3) insufficient sampling for accurate derivative evalua- 
tions. Nodes are usefully divided into two categories;^ 
depending on whether Q is formally well-behaved ( "type 
one" nodes) or singular ("type two" nodes). For sta- 
tionary state solutions to the Schrodinger equation, for 
instance, all nodes are type one nodes. In principle, type 
one nodes are "gentler" than type two nodes; however, 
from a numerical standpoint, even type one nodes will 
give rise to the problems listed above, because the slight- 
est numerical error in the evaluation of Q is sufficient to 
cause instability. 

In the best case, the node problem simply results in 
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substantially more trajectories and time steps than the 
corresponding classical calculation; in the worst case, the 
QTM calculation may fail altogether, beyond a certain 
point in time. Several numerical methods, both "exact" 
and approximate, are currently being developed to deal 
with this important problem. The latter category in- 
cludes the artificial viscosity'^^-'^^ and linearized quantum 
force methods, both of which have proven to be very 
stable. While such approximate methods may not cap- 
ture the hydrodynamic fields with complete accuracy in 
nodal regions, they do allow for continued evolution and 
long-time solutions, often unattainable via use of a tradi- 
tional QTM. The "exact" methods include the adaptive 
hybrid methods,— and the complex amplitude method.— 
In the adaptive hybrid methods, for which hydrodynamic 
trajectories are evolved everywhere except for in nodal 
regions, where the time-dependent Schrodinger equation 
is solved instead to avoid node problems. Although they 
have been applied successfully for some problems, these 
methods are difficult to implement numerically, since not 
only must the hydrodynamic fields be somehow moni- 
tored for forming singularities, but there must also be 
an accurate means for interfacing and coupling the two 
completely different equations of motion. The complex 
amplitude method is cleaner to implement, but is only 
exact for linear and quadratic Hamiltonians. 

In a recent paper^^i hereinafter referred to as "pa- 
per I," one of the authors (Poirier) introduced a new 
strategy for dealing with the node problem, based on a 
bipolar decomposition of the wavefunction. The idea is 
to partition the wavefunction into two (or in principle, 
more) component functions, i.e. ^E" = -I- ^"2. One 
then applies QTM propagation separately to and '^2, 
which can be linearly superposed to generate 5" itself at 
any desired later time. In essence, this works because the 
Schrodinger equation itself is linear, but the equivalent 
Bohmian mechanical, or quantum Hamilton's equations 
of motion (QHEM) are not.^^ In principle, therefore, one 
may improve the numerical performance of QTM cal- 
culations simply by judiciously dividing up the initial 
wavepacket into pieces. 

Although bipolar decompositions have been around 
for quite some time^i^ their use as a tool for cir- 
cumventing the node problem for QTM calculations 
is quite recent. Two promising new exact methods 
that seek to accomplish this are the so-called "counter- 
propagating wave" method (CPWM),^^ and the "cov- 
ering function" method (CFM).4? In the CPWM, the 
bipolar decomposition is chosen to correspond to the 
semiclassical WKB approximation,^"* for which all of the 
hydrodynamic field functions are smooth and classical- 
like, and the component wavefunctions are node-free. 
Interference is achieved naturally, via the superposition 
of left- and right-traveling (i.e. positive- and negative- 
momentum) waves. For one-dimensional (ID) stationary 
bound states, it can be shown that the resultant bipo- 
lar quantum potential q{x) becomes arbitrarily small in 
the large action limit, even though the number of nodes 



becomes arbitrarily large. (Note: in accord with the con- 
vention established in Ref. [U, upper/lower case will be 
used to denote the unipolar/bipolar field quantities). In 
the CFM, the idea is to superpose some well-behaved 
large-amplitude wave, with the actual ill-behaved (nodal 
or wildly oscillatory) wave, so as to "dilute" the undesir- 
able numerical ramifications of the latter. 

This paper is the second in a series designed to ex- 
plore the CPWM approach, introduced in paper I. As 
discussed there in greater detail, there are many motiva- 
tions for this approach, but the primary one is to recon- 
cile the semiclassical and Bohmian theories, in a manner 
that preserves the best features of both, and also sat- 
isfies the correspondence principle. For our purposes, 
this means that the Lagrangian manifolds (LMs) for the 
two theories should become identical in the large action 
limit (Sec. Ill Bp . As described above, a key benefit of 
the CPWM decomposition is an elegant treatment of in- 
terference, the chief source of nodes and "quasi-nodes" — 
(i.e. rapid oscillations) in quantum mechanical systems. 
An interesting perspective on the role of interference in 
semiclassical and Bohmian contexts is to be found in a 
recent article by Zhao and Makrii^ 

Whereas paper I focused on stationary bound states for 
ID systems, the present paper (paper II) and the next 
in the series (paper III)i^ concern themselves with sta- 
tionary scattering states. The CPWM decomposition of 
paper I is uniquely specified for any arbitrary ID state — 
bound or scattering — and in the bound case, always sat- 
isfies the correspondence principle. However, the non-L^ 
nature of the scattering states is such that the paper I 
decomposition generally does not satisfy the correspon- 
dence principle in this case. Simply put, the quantum 
trajectories and LMs exhibit oscillatory behavior in at 
least one asymptotic region (thereby manifesting reflec- 
tion), whereas the semiclassical LMs do not. This is not 
a limitation of the CPWM, but is rather due to the fun- 
damental failure of the basic WKB approximation to pre- 
dict any reflection whatsoever for above-barrier energies, 
as has been previously well established^ii^i^ In semi- 
classical theory, a modification must therefore be made to 
the basic WKB approximation, in order to obtain mean- 
ingful scattering quantities. As discussed in Sec. Ill Bl and 
in paper HI, our approach will be to apply a similar mod- 
ification to the exact quantum decomposition (actually, a 
reverse modification) such that the correspondence prin- 
ciple remains satisfied, and the two theories thus recon- 
ciled, even for scattering systems. 

It will be shown the modified CPWM gives rise to 
bipolar Bohmian LMs that are identical to the semi- 
classical LMs, regardless of whether or not the action 
is large. Put another way, this means that the bipolar 
quantum potentials q effectively vanish, so that the resul- 
tant quantum trajectory evolution is completely classical. 
Moreover, the resultant component wavefunctions, ^'i(x) 
and (x), correspond asymptotically to the familiar "in- 
cident," "transmitted," and "reflected" waves of tradi- 
tional scattering theory. Thus, the modified CPWM im- 
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plementation of the bipolar Bohmian approach provides 
a natural generalization of these conceptually fundamen- 
tal entities throughout all of configuration space, not just 
in the asymptotic regions, as is the case in conventional 
quantum scattering theory. 

The above conclusions will be demonstrated for both 
discontinuous and continuous potential systems, in pa- 
pers II and III, respectively. Discontinuous potentials — 
e.g. the hard wall, the step potential, and the square 
barrier/well — serve as a useful benchmark for the modi- 
fied CPWM approach, because the scattering component 
waves (e.g. "incident wave," etc.) in this case are well- 
defined throughout all of configuration space, according 
to a conventional scattering treatment. Although this is 
no longer true for continuous potentials, the foundation 
laid here in paper II can be extended to the continu- 
ous (and also time-dependent) case as well, as described 
in paper III. Additional motivation for the development 
of a scattering version of the CPWM, vis-a-vis the rel- 
evance for chemical physics applications, is provided in 
paper III. Additional motivation for the consideration of 
discontinuous potentials is provided in Sec. Ill Bl of the 
present paper. 



II. THEORY 



equations, i.e. 



A. Background 

1. Bohmian mechanics 



(4) 



the continuity equation ensures that the probability [i.e. 
density, R{x,t)^, times volume element] carried by indi- 
vidual quantum trajectories is conserved over the course 
of their time evolution. 



2. CPWM decomposition for stationary states 



In paper I, we derived a unique bipolar decomposition. 



(5) 



for stationary eigenstates '^{x) of ID Hamiltonians of the 
Eq. ^ form, such that: 

1. 5'±(x) are themselves (non-L^) solutions to the 
Schrodinger equation, with the same eigenvalue, E, 
as ^'(x) itself. 

2. The invariant flux values, ±_F, of the two solutions, 
\E'±(a;), equal those of the two semiclassical (WKB) 
solutions. 

3. The median of the enclosed action, xq, equals that 
of the semiclassical solutions. 



According to the Bohmian formulation,— the 
QHEM are derived via substitution of the ID (unipolar) 
wavefunction ansatz, 



«'(x,t) = i?(a;, i)e*^(^'*'/'^ 



(1) 



into the time-dependent Schrodinger equation. For the 
ID Hamiltonian, 



2m dx^ 



(2) 



this results in the coupled pair of nonlinear partial dif- 
ferential equations. 



dS{x,t) 

dt 
dR{x, t) 

dt 



R" 



2m ^(^)+2mi? 

= ^Lr' s' -^RS", 
m 2m 



(3) 



where m is the mass, V{x) is the system potential, and 
primes denote spatial partial differentiation. 

The first of the two equations above is the quantum 
Hamilton- Jacobi equation (QHJE), whose last term is 
equal to —Q{x,t), i.e. comprises the quantum potential 
correction. The second equation is a continuity equation. 
When combined with the quantum trajectory evolution 



There are other important properties of the ^±{x)~ as 
discussed in Sec. U and in Ref. '2^. Nevertheless, the 
above three conditions are sufficient to uniquely specify 
the decomposition. In the special case of bound (i.e. L^) 
stationary states, the real-valuedness of 5* (a;) implies that 
the '^±{x) are complex conjugates of each other. 



B. Scattering systems 

It is natural to ask to what extent the above analysis 
may be generalized for scattering potentials. Certainly, 
^(x) itself is no longer L^, nor even real- valued, and 
there are generally two linearly independent solutions 
of interest for each E, instead of just one. Condition 
(1) above poses no difficulty for 5'±(x), as these com- 
ponent functions are non-L^ and complex-valued, even 
in the bound eigenstate case. In principle, condition (2) 
is not difficult either; although the flux value depends 
on the normalization of 5* itself, which is not L^, cer- 
tain well-established normalization conventions for scat- 
tering states exist, that can be applied equally well to 
semiclassical and exact quantum solutions. There is no 
action median per se for scattering states, as the action 
enclosed within the '^+{x) and ^'_(a;) phase space La- 
grangian manifolds— (LMs) is infinite; however, 
the scattering analog of condition (3) is related to the 



4 



asymptotic boundary conditions, and it is here that one 
encounters difficulty. Moreover, an additional concern is 
raised by the doubly-degenerate nature of the continuum 
eigenstates, namely: should each scattering ^(a;) have 
its own 4'±(x) decomposition, or should there be a sin- 
gle ^'±(a;) pair, from which all degenerate ^'(x)'s may be 
constructed via arbitrary linear superposition? 

To resolve these issues, we will adopt the same gen- 
eral strategy used in paper I, i.e. we will resort to semi- 
classical theory as our guide, wherever possible. We 
will also exploit certain special features of the scatter- 
ing problem not found in generic bound state systems, 
such as the asymptotic potential condition V'{x) —^ 
as a; — > ±00 (where primes denote spatial differentia- 
tion), and its usual implications for scattering theory and 
applications.^"'^ 

The basic WKB solutions are given by 

vI/^±^(a;)=r,e(a;)e±^^-(^)/^ (6) 

where 



rscix) = J ^-^ and = V2m [E ~ V{x)] 

(7) 

The corresponding positive and negative momentum 
functions, specifying the semiclassical LMs, are given by 
p±{x) = ±s'g^{x). Equations © and ([7]) apply to both 
bound and scattering cases; note that for both, ^Pj? (x) are 
complex conjugates of each other. The asymptotic poten- 
tial condition ensures that these approach exact quantum 
plane waves asymptotically, with the usual scattering in- 
terpretations, i.e. in the x — > —00 asymptotic 
region is the incoming wave from the left (usually taken 
to be the incident wave), ^'+(x) as x ^ 00 is the out- 
going wave from the left (the usual transmitted wave), 
etc. 

Insofar as determining the corresponding exact quan- 
tum solutions '9±{x), the procedure described in paper 
I is still appropriate for bound and semi-bound (i.e. on 
one side only) states, in that the results satisfy the corre- 
spondence principle globally, as desired (for semi-bound 
examples, consult the Appendix). For true scattering 
states, however, this procedure fails, in the sense that if 
^+{x) is chosen to match the normalization and flux of 
^^(x) in the x 00 asymptote, then it will necessar- 
ily approach a nontrivial linear superposition of ^'^(a;) 
and 4'!!:(a;) in the x — *■ —00 asymptote, and vice- versa. 
There is therefore an ambiguity as to how the corre- 
sponding quantum \E'±(a;)'s should be defined, i.e. which 
asymptotic region should be used to effect the correspon- 
dence. More significantly though, either choice will result 
in component functions '^±{x) with substantial interfer- 
ence in one of the two asymptotic regions. This is due to 
partial reflection of the exact quantum scattering states, 
which is not predicted by the basic WKB approxima- 
tion. Thus, in the large action limit, the exact quan- 
tum solutions manifest large-magnitude quantum poten- 
tials, q±{x), and rapidly oscillating field functions q±{x), 



r±{x)^ and p±{x) — exactly the undesirable behavior that 
the CPWM was introduced to avoid — whereas the corre- 
sponding basic WKB functions are smooth, and asymp- 
totically uniform. 

The lack of any partial reflection is a well-understood 
shortcoming of the WKB approximation^^i^i^i^^ — i.e., 
the basic ^'!f(a:) components, though elegantly con- 
structed from smooth classical functions rsc{x) and 
Ssc{x), do not in and of themselves correspond to any 
actual quantum scattering solutions ^>{x). In light of the 
bipolar decomposition ideas introduced in paper I, how- 
ever, our perspective is the reverse one: for any actual 
quantum ^'(x), can one determine an Eq. ^ decompo- 
sition such that the resultant ^±{x) resemble their well- 
behaved semiclassical counterparts, and is such a decom- 
position unique? Among other properties,— the Vl'±(x) 
LM's should become identical to the semiclassical LM's in 
the large action limit, so as to satisfy the correspondence 
principle. Based on the considerations of the previous 
paragraph it is clear that the paper I decomposition does 
not achieve this goal, when applied to stationary scatter- 
ing states. 

We defer a full accounting of these issues — in the 
context of completely arbitrary continuous potentials 
V{x) — to paper III, wherein it will be demonstrated how 
to compute exact quantum reflection and transmission 
probabilities (and stationary scattering states) using only 
classical trajectories, and without the need for explicit 
numerical differentiation of the wavefunction. In the 
present paper, we lay the foundation for paper III, by 
focusing attention onto two key aspects whose develop- 
ment comprises an essential prerequisite. 

First, as the paper III approach treats V{x) as a se- 
quence of stepS)^ the present paper II will focus exclu- 
sively on the step potential and related discontinuous 
potential systems, for which V(x) = const in between 
successive steps. Discontinuous potentials are important 
for chemical physics, because they model steep repulsive 
wells, and are used in statistical theories of liquids. More- 
over, they hold a special significance for QTM methods, 
for which they serve as a "worst-case scenario" bench- 
mark. Indeed, conventional QTM techniques always fail 
when applied to discontinuous potentials. To date. The 
only such calculations that have been performed^ have 
computed the quantum potential from a completely sep- 
arate time-dependent fixed-grid calculation (the "analyt- 
ical approach")!* rather than directly from the quan- 
tum trajectories themselves. Even if one could propagate 
trajectories for discontinous systems using a traditional 
QTM, the trajectories that would be generated would be 
very kinky and erratic)^ and a great many time trajec- 
tories and time steps would thus be required. 

Second, since the new ^±{x) do not satisfy condition 
(1), unlike the paper I CPWM decomposition, the time 
evolution of these two component functions is clearly not 
that of the time-dependent Schrodinger equation. More- 
over, since the |^'±P are constant over time [because 
itself is stationary, and Eq. ([5]) is presumed unique], the 
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two "^±{x,t) time evolutions must be coupled together. 
It is essential that the nature of this coupling be com- 
pletely understood, in order that the present approach 
may be generalized to non-stationary state situations — 
e.g. to wavepacket scattering, as will be discussed in 
future publications. The ramifications for QTMs are 
equally important. Accordingly, the present paper fo- 
cuses on the QTM propagation of the wavefunction and 
its bipolar components — with a keen eye towards gener- 
ality and physical interpretation — even though the states 
involved are stationary. This approach leads to a peda- 
gogically useful reinterpretation of "incident," "transmit- 
ted," and "reflected" waves — very reminiscent of ray op- 
tics in electromagnetic theory — which is applicable much 
more generally than traditional usage might suggest. 

C. Basic applications 

The necessary theory will be developed over the course 
of a consideration of various model application systems 
of increasing complexity. 

1. free particle system 

Let us first consider the simplest case imaginable, the 
free particle system, V{x) = 0. In this case, the exact so- 
lutions ^'±(x) = ^!f(a:) clearly satisfy the conditions of 
Sec. Ill A"2l and the bipolar quantum potentials q±{x) are 
zero everywhere. Thus, the bipolar decomposition devel- 
oped for bound states in paper I can be used directly with 
this continuum system, requiring only the slight modifi- 
cation that arbitrary linear combinations of ^'^(a;) and 
^'f5(a;) are to be allowed, in order to construct arbitrary 
scattering solutions '^(x). For convenience, the linear 
combination coefficients will from here on out be directly 
incorporated into the amplitude functions, r±{x), and 
phase functions, s±{x), so that Eq. (O is still correct. 

If from all solutions ^'(x) one considers only that which 
satisfies the usual scattering boundary conditions (i.e. in- 
cident wave incoming from the left) then the negative mo- 
mentum wave vanishes, and ^!{x) = There 
is zero reflection, and 100% transmission. Put another 
way, the incident flux, limx^-oc is equal to the 

transmitted flux, Imix^+Qc j+{x), where 

■ / N ^ Tt*/ ^d^±{x) d^*^{x) ' 
Aim [ dx dx 

[both flux values are equal to F, as in Eq. ([7])]. 

In the quantum trajectory description, flux mani- 
fests as probability-transporting trajectories, which move 
along the LMs. For the boundary conditions described 
above, there are only positive momentum trajectories, 
moving uniformly from left to right with momentum 



p+{x) = \/2mE. If a -{x) contribution were present, 
its trajectories would move uniformly in the opposite di- 
rection [p^{x) = —\/2mE] Since the two components 
±{x) are in this case uncoupled, the positive and neg- 
ative momentum trajectories would have no interaction 
with each other. 



2. hard wall system 
We next consider the hard wall system: 

Vix) = I ^ for X < 0; 

^^^> loo forx>0. ^' 

In the X < region, the two ^±{x) components are ex- 
actly the same as in the free particle case, except that 
the ^'(0) = boundary condition imposes the additional 
constraints, 

s_(0) = s+(0) + 7rmod(27r) ; r_(0)=r+(0). 

(10) 

This also results in only one linearly independent solution 
instead of two, i.e. ^'(x) oc sin(fcx), with k = ^/2mE/'h. 
Regarding the LMs and trajectories, in the x < region, 
these are identical to those of Sec. Ill CTl e.g. the ^+(x) 
LM trajectories move uniformly to the right, towards the 
hard wall at x = 0. 

It is natural to ask what happens when the ^'+(x) LM 
trajectories actually reach x = 0. There are two rea- 
sonable interpretations. The first is that the trajectories 
keep moving uniformly into the x > region of configura- 
tion space. This approach treats the hard wall system as 
if it were the free particle system, but with the x > re- 
gion effectively ignored;^ This underscores the fact that 
unlike ^'(x) itself, the individual 5'±(x) components per 
se are unconstrained at the origin — though the Eq. (jlOp 
constraint implies a unique correspondence between the 
two. This interpretation also makes it clear that for the 
hard wall system, the paper I decomposition is essentially 
identical to the present decomposition, as is worked out 
in detail in the Appendix. 

In the second interpretation, the effect of the hard 
wall at X = is to cause instantaneous elastic reflec- 
tion of a \E'+(x) LM trajectory momentum, from p = 
p+ = +\/2mE to p — p^ = —\/2mE. Afterwards, 
the reflected trajectory propagates uniformly backward, 
along the \l/_(x) LM. In this interpretation, the trajec- 
tories never leave the allowed configuration space, x < 0. 
However, wavepacket reflection is essentially achieved via 
trajectory hopping from one LM to the other — not un- 
like that previously considered, e.g., in the context of 
non-adiabatic transitions.-54 The trajectory hopping in- 
terpretation is adopted in the present paper, and in pa- 
per III, but the first interpretation will also be recon- 
sidered in later publications. Note that for discontin- 
uous potentials — and indeed more generally^^ — one can 
regard trajectory hopping as the source of ^'±(x) inter- 
action coupling. 
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For the hard wall case, trajectory hopping only mani- 
fests at X = 0, the sink of all '^+{x) LM trajectories, and 
the source of all LM trajectories. If these trajec- 

tories are to be regarded as one and the same via hop- 
ping, then a unique field transformation for r, s, and all 
spatial derivatives, must be specified. Fortunately, the 
unique correspondence between and ^'_(a;) de- 

scribed above, enables one to do just that. In particular, 
Eq. pO)) specifies the correct transformations for r and s, 
as transported by the quantum trajectories. All spatial 
derivatives of arbitrary orders can then be obtained via 
spatial differentiation of Eq. ^ — although in the hard 
wall case, only the s' condition, p-{0) = —p+{0) is rele- 
vant, because all higher order derivatives are identically 
zero. 

Since the magnitudes of the p and r fields associated 
with a given quantum trajectory are unchanged as a re- 
sult of the trajectory hop, Eq. ([8]) implies that the inci- 
dent and reflected flux values are the same (apart from 
sign), and so the scattering system exhibits 100% reflec- 
tion and zero transmission (along each LM, the flux is 
invariant^^) . These basic facts of the hard wall system are 
of course well understood. The point, though, is that we 
have now obtained the information in a time- dependent 
quantum trajectory manner, rather than through the 
usual route of applying boundary conditions to time- 
independent piecewise component functions. In other 
words, Eq. (|10|) now refers to individual quantum tra- 
jectories, rather than to wavefunctions. 

This shift of emphasis is very important, and leads 
to quite a number of conceptual and computational ad- 
vantages. For instance, the standard description of the 
hard wall stationary states would decompose these into 
plane wave components interpreted as "incident" and 
"reflected" waves. This language suggests a process, or 
change over time — i.e. a state that is initially incident, 
at some later time is somehow transformed into a re- 
flected state. Nothing in the standard description, how- 
ever, would seem to render transparent the usage of such 
terminology, i.e. ^!{x) is stationary, and so the reflected 
and transmitted components are in fact both present for 
all times. Of course, a localized superposition of station- 
ary states, i.e. a wavepacket, may well exhibit such an 
explicit transformation over the course of the time evo- 
lution, as such a state is decidedly non-stationary. In- 
deed, wavepackets are relied upon by the more rigorous 
formulations of scattering theory, in order to justify the 
use of terms such as "reflected wave," even in a station- 
ary context.'^ Such formulations, though certainly legit- 
imate, seem always to require a clever use of limits, the 
subtle distinction between unitary and isometric trans- 
formations, and other esoteric mathematical tricks. 

On the other hand, the time-dependent bipolar quan- 
tum trajectory hopping picture presented above provides 
a physicality to such language that is immediately ap- 
parent. Over the course of the time evolution, although 
the wavefunction as a whole is stationary, each individ- 
ual trajectory is first incident from the left, then collides 



with the hard wall, and is subsequently reflected back 
towards the left (i.e. towards x —^ —oo). The bipo- 
lar quantum trajectories are all classical, as the bipolar 
quantum potentials, q±{x), are zero everywhere except 
at the wall itself. Interference arises naturally from the 
superposition of the two LMs — i.e., from the trajectories 
that have already progressed to the point of reflecting, 
vs. those that have not reflected yet. In contrast, since 
^'(x) itself exhibits very substantial interference, and an 
infinite number of nodes, the traditional unipolar QTM 
treatment would be very ill-behaved, i.e. R{x) would 
oscillate wildly in the large k limit, and Q{x) would be 
numerically unstable near the nodes. Apart from these 
important pragmatic drawbacks, the incident/reflected 
interpretation of the quantum trajectories would also be 
lost. 

The bipolar quantum trajectory description of the hard 
wall system is very reminiscent of ray optics, as used to 
describe the reflection of electromagnetic waves off of a 
perfectly reflecting surfacci^ Indeed, much can be gained 
from applying a ray optics analogy to quantum scattering 
applications, especially where discontinuous potentials 
are concerned. One can construct a simple gedanken- 
experiment as follows. Let xl < denote some effective 
left edge of the system, well to the left of the interaction 
region. At some initial time t = 0, all trajectories on the 
positive LM lying to the right of xl are ignored, as is the 
negative LM altogether. One then evolves the retained 
trajectories over time, and monitors the contribution that 
just these trajectories make to the total wavefunction. In 
some respects, it is as if the point xl were serving as the 
initial wavefront for some incoming wave, that at t = 
had not yet reached the hard wall/reflecting surface. Of 
course, if the actual wave were in fact truncated in this 
fashion, then the discontinuity in the field functions at 
the wavefront would result in a very non-trivial propaga- 
tion over time, owing to the high-frequency components 
implicitly present. For actual waves, the precise nature 
of the wavefront is known to have a tremendous impact 
on the resultant dynamics.— We avoid such complicat- 
ing details by always interpreting the "actual wave" to 
be the full stationary wave itself, i.e. the truncation is 
conceptual only. 

In the ray optics analogy, the above situation is like a 
source of light located at , which is suddenly "turned 
on" at i = 0. It takes time for the wavefront to propagate 
to the reflecting surface, and additional time for the re- 
flected wavefront to make its way back to x — x^. Prior 
to the latter point in time, the evolution of the trun- 
cated electromagnetic wave is decidedly non-stationary; 
afterwards however, a stationary wave is achieved, at 
least within the region of interest, xl ^ x < 0, as the 
wavefront has by this stage propagated beyond this re- 
gion. The same qualitative comments apply to the bipo- 
lar quantum case, although of course the evolution equa- 
tions are different. 

A similar prescription may be used to achieve rudi- 
mentary "wavepacket dynamics," even in the context of 
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purely stationary states. Instead of retaining all initial 
trajectories that lie to the left of a;^, one retains only 
those that lie within some finite interval. The result- 
ing time evolution is analogous to a light source that is 
turned on at i = 0, and then turned off at some later time 
(prior to when the wavefront arrives at the reflecting sur- 
face) . The initial "wavepacket" has uniform density, and 
moves with uniform speed towards the hard wall. Inter- 
ference fringes then form after the foremost trajectories 
have been reflected onto the negative LM. Eventually, 
all trajectories within the interval are reflected, at which 
point interference ceases (the nodes are "healed"—), uni- 
form density is restored, and the reflected wave travels 
with uniform speed in the reverse direction, back towards 
the starting point xl- Qualitatively, this behavior is 
clearly similar to that undergone by actual wavepackets 
reflecting off of barrier potentials. 



D. More complicated applications 

The ideas described above can be easily extended to 
more complicated discontinuous potential systems, such 
as up- and down-step potentials, and any combination of 
multiple steps, e.g. square barriers and square wells. In 
paper III, they will even be extended to arbitrary contin- 
uous potentials.— In every case, the ray optics analogy 
from electromagnetic theory may also be extended ac- 
cordingly. This approach provides a useful perspective on 
global reflection and transmission in scattering systems, 
and in particular, demonstrates how such quantities may 
be obtained from a single, universal expression for local 
reflection and transmission. 



1. step potential system — above barrier energies 



We next consider the step potential system: 



Vix) 





Vo 



for a; < 0; 
for X > 0, 



(11) 



Classically, this system exhibits 100% transmission if the 
trajectory energy is above the barrier (i.e. E > Vq), and 
100% reflection if the trajectory energy is below the bar- 
rier {E < Vo)- Quantum mechanically, all above barrier 
trajectories are found to exhibit partial reflection and 
partial transmission, although there is a general increase 
in transmission probability with increasing energy. The 
below barrier quantum trajectories exhibit 100% reflec- 
tion, as in the classical case; however, they also mani- 
fest tunneling into the classically forbidden a; > region. 
Thus even quantum mechanically, the the above and be- 
low barrier cases must be handled somewhat differently. 

To begin with, we consider the above-barrier case. 
Note that the LM's are unbounded in either direction, i.e. 
the classically allowed region extends to both asymptotes, 
X ±oo. Incoming trajectories can therefore originate 



from either asymptote, thus giving rise to two linearly 
independent solutions, ^'(a;). This is in stark contrast 
to the hard wall system, for which incoming trajectories 
could only originate from x — *■ — oo, thus resulting in only 
one linearly independent solution for 5" (a;). 

In the standard time-independent picture, one starts 
with the four piecewise solutions. 



and *B±(a;) = e^^P^^/'', (12) 



where region A corresponds to a; < 0, region _B to a: > 0. 
The momenta values are classical, i.e. 



PA 



V2mE and pB = y^2m{E-Vo). (13) 



Matching '^{x) and ^'(x) boundary conditions at a; = 0, 
and specifying asymptotic boundary conditions for ^'(a;), 
then enables a unique determination of the four complex 
coeflicients A± and B± in 

M/M - /"^+*^+(^)+^-*^-(^) fora;<0; , . 
\B+*b+(.t) +B_*B-(a;) for a; > 0, ' ^ ' 

In general, the solution coefficients depend on the par- 
ticular stationary solution of interest. For the usual scat- 
tering convention of an incident wave incoming from the 
left (Fig. [T|) the solutions are 



B.=T = 



A+ = l 

2pA 



PA + PB 



A^ =R = 
S_ = 0, 



/ PA~PB 
\PA +PB 



(15) 



where R and T are respectively, refiection and transmis- 
sion amplitudes. When fiux is properly accounted for, 
the resultant reflection and transmission probababilities 
(which add up to unity) are given by 



m 



Pt, 



PA 



(16) 



Note that Eqs. (fT5|) and p6|) above are correct for both 
an "up-step" and a "down-step" — i.e. for Vq positive or 
negative. We can also apply these equations to the "op- 
posite" boundary conditions, i.e. to an incident wave 
incoming from the right, by simply transposing A and B 
subscripts, and + and — subscripts {pA and pb are still 
positive). This is important, because any stationary so- 
lution 4" (a;) can be obtained as some linear superposition 
of left-incident and right-incident solutions. 

Regarding the time-dependent interpretation, it is ev- 
ident that upon reaching the step discontinuity, left- 
incident trajectories must be partially reflected and par- 
tially transmitted. The trajectory is suddenly split into 
two, one that continues to propagate along the positive 
LM for the transmitted B region (i.e. the B+ LM) and 
the other being instantaneously reflected down to the 
A— LM. Moreover, since probability carried by individ- 
ual quantum trajectories is conserved, ^^'"^^ this splitting 
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X (a.u.) 

FIG. 1: Component waves for a left-incident stationary eigen- 
state of the up-step barrier problem with £ > Vb, as described 
in Sec. IIV Al Solid and dashed lines represent real and imagi- 
nary contributions, respectively. The dot-dashed line denotes 
the location of the step. 

must be done in a manner that preserves both probability 
and flux. In other words, the local splitting of the trajec- 
tory at a; = must correspond to Eq. (jlSp . which is now 
regarded as a local condition, giving rise to local reflection 
and transmission amplitudes, R and T. For the present 
step potential case, these local quantities are directly re- 
lated to the global P^cd and Ptrans values via Eq. 
For multiple step potentials (Sec. HID 3p . the global ex- 
pressions above [Eq. p6)) ] no longer apply; however, a 
local, time-dependent trajectory version of Eq. psp does 
turn out to be correct. 

Such an expression, immediately applicable to all sin- 
gle and multiple step systems, can be written as follows: 

n-cfi = 1 ^'Z*^ Ptrani, \ ^^^^^ ^ ^^^^ ^ ^^^^^ ^^^^ 
VPi/r + Ptrans/ 

^trans — I , I ^inc ; Strans — 5inc-l,-Loj 

\Pi/r + Ptrans/ 

In the above equations, "inc" refers to any trajectory, 
locally incident on some particular step from some par- 
ticular direction, which spawns both a locally reflected 
trajectory, "refl," and a locally transmitted trajectory, 
"trans" . The quantity pi/^ is the (positive) momentum 
associated with the locally incident/reflected trajectory; 
similarly, ptrans (also positive) is associated with the lo- 
cally transmitted trajectory. For above-barrier incident 
trajectories, note that the local reflection and transmis- 
sion amplitudes are both real, thus ensuring the reality 
of r and s for the spawned trajectories. 

Returning to the step potential system, the ray optics 
picture can once again shed some interesting light. The 
optical analog of the step is an interface between two 



media with different indices of refraction. Light incident 
on such an interface will partially reflect back towards 
the original source, and partially refract forwards into 
the new medium. The refraction is completely analogous 
to the discontinuous change in momentum, {ps — Pa), 
that suddenly occurs as one crosses the step (Fig. [3]). In 
any event, the ray optics gedankenexperiment described 
in Sec. Ill C 21 can also be applied to the step potential 
system, in order to obtain a particular stationary solution 
"^{x) with any desired boundary conditions. 

For instance, suppose one is interesting in constructing 
the left-incident wave solution, i.e. that of Eq. (fT5| . At 
t = 0, only the '^a+ wave is considered, and only those 
trajectories for which x < xl, as before. As the incident 
trajectories reach the step, two new waves are dynami- 
cally created from the spawned trajectories: a transmit- 
ted wave traveling to the right, and a reflected wave trav- 
eling to the left. A plot of the overall density |^'(a;, t)p so 
obtained will change over time, as the transmitted and re- 
flected wavefronts propagate into their respective regions 
(Fig. [5]). Eventually, however, these wavefronts will prop- 
agate beyond the region of interest, i.e. xl < x < xr, 
where > is the right edge of the region of inter- 
est. When this occurs, the solution for \E'(x) obtained 
within the region of interest will be exactly equal to the 
stationary solution with the desired boundary condition. 

As in the hard wall case, one can also perform step 
potential "wavepacket dynamics" by restricting consid- 
eration to just those initial a+ trajectories lying within 
some coordinate interval. The wavepacket will propagate 
towards the step with uniform density and speed. As the 
flrst few trajectories hit the step, a uniform transmitted 
wave will be formed in the B region. In the A region, 
the sudden appearance of a a- wave will introduce in- 
terference wiggles in the overall density plot (although 
no nodes per se, owing to partial reflection only). Even- 
tually, after all trajectories have progressed beyond the 
step, well-separated transmitted and reflected wavepack- 
ets emerge, propagating in their respective spaces and di- 
rections. There is no longer any interference in the A re- 
gion, as the incident wave is now gone, having been com- 
pletely divided into the two final contributions. Pj-cfl and 
Ptrans valucs may be determined via monitors placed at 
XL and Xr, either by integrating probability over time as 
the respective wavepackets travel through, or by record- 
ing the (constant) amplitude values R and T, and apply- 
ing Eq. p6|) . 



2. step potential system — below barrier energies 

The case for which the incident trajectory energies are 
below Vq requires special discussion. In this case, the 
classical LMs and trajectories are confined to the A re- 
gion only (i.e. to x < 0), as the entire B region is classi- 
cally forbidden. In the language of Sec. Ill Bl these below 
barrier states are therefore semi-bound, implying that 
there is only one linearly independent stationary solu- 
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X (a.u.) 

FIG. 2; Wavefunction plot for a left-incident stationary eigen- 
state of the up-step barrier problem with E <Vo, as described 
in Sec. IIV Al Solid and dashed lines represent real and imag- 
inary contributions, respectively, for the analytical solution. 
Squares and circles denote corresponding numerical results. 
The shaded box represents the tunneling region. 



tion, ^'(x) which without loss of generality, must be real- 
valued. This in turn implies that the 5'±(x) are complex 
conjugates of each other, as in the bound state case dis- 
cussed in paper I. Indeed, one option is to simply apply 
the paper I decomposition to such problems. This ap- 
proach is discussed in detail in the Appendix, wherein it 
is shown to provide a natural extension of classical tra- 
jectories into the tunneling region. 

On the other hand, the trajectory hopping-based de- 
composition scheme offers a different, but also very nat- 
ural means to accomplish the same task — which has the 
added advantage that all bipolar quantum potentials van- 
ish, except at a: = 0. The idea is simply to treat all ex- 
pressions in Sec. HID II as being literally correct for the 
below barrier case as well, with the understanding that 
the requisite quantities need no longer be real. In partic- 
ular, ptrans — [Eq. (HH)] bccomcs purc positive imag- 
inary, implying that the transmitted trajectories "turn a 
corner" in the complex plane, and start heading off in the 
positive imaginary direction, with speed tiK/m (Fig. [S]). 
Along this path, the transmitted wave is an ordinary 
plane wave; however, when analytically continued to the 
real axis in the a; > region (via a 90° clockwise rotation 
in the complex plane), the familiar exponentially damped 
form results (Fig. ^ . 

For the reflected wave, Eq. p?|) states that the re- 
flected "phase" remains unchanged. However, rj-cn is 
now complex, leading to an effective phase shift of 2d, 
where 6 is defined in the Appendix [Eq. (HI])]. For lo- 
calized wavepackets, a physical significance can be at- 
tributed to this phase shift, in both quantum mechanics 
and electromagnetic theory; it is the source of the Goos- 
Hanchen effectf^!^ a time delay observed in conjunction 



with total internal reflection. Consequently, in the time- 
dependent wavepacket context, it may be more appro- 
priate to associate the phase shift with time, rather than 
with s or r — specifically, with the delay time needed to 
accrue sufficient action so as to compensate for the shift. 
For stationary states, however, such a time delay would 
be inconsequential, because all trajectories are identical 
apart from overall phase. Consequently, we do not con- 
sider such time delays explicitly in this paper, though we 
will return to this issue in future publications. 



3. multiple step systems 

The most interesting case is that for which there are 
multiple discontinuities, occurring at arbitrary locations 
Xk (with k = 1,2, ... ,1), and dividing up configuration 
space into / + 1 regions, labeled A, B, C , etc. In each re- 
gion, the potential energy has a different constant value, 
i.e. V{x) = Va in region A, etc. From an optics point 
of view, this system is analogous to a stack of different 
materials, each with its own thickness, and index of re- 
fraction. Our primary focus in this paper will be square 
barrier/well systems for which 1 = 2, and Va = Vc. How- 
ever, all of the present analysis extends to the more gen- 
eral case described above. 

In the standard time-independent picture, the solu- 
tion is obtained via a straightforward generalization of 
Eqs. (Uni), and p^ . However, even when compara- 

ble left-incident boundary conditions are specified as in 
Sec. iHDTI -i.e. A+ = l', and (for I = 2) C- = 0— the 
remaining coefficient values are fundamentally different 
from those of the single-step case. To begin with, only 
the rth step exhibits the characteristics of a (locally) left- 
incident single-step solution; all other steps involve four 
non-zero coefficients, corresponding locally to some su- 
perposition of left- and right-incident waves. Even more 
importantly, however, the expressions for the coefficient 
values as a function of system parameters in no way re- 
sembles Eq. (dSI; in particular, these now depend explic- 
itly on the Xk values, as well as on Va, Vb, etc. The same 
is also true for the global Pj-ofl and Ptrans expressions, as 
compared with Eq. p^ . 

It is this dependence on the other steps that gives rise 
to the global nature of the time-independent solutions; 
i.e. the coefficient values at one step depend in principle 
on the properties of all of the other steps, no matter how 
far away these might be located. Consequently, a reflec- 
tion probability as obtained from the A^ value associated 
with the first, fc = 1 step, cannot be determined without 
extending the analysis out to the final step at x = xi, 
in the standard time-independent picture. On the other 
hand, a primary goal of the time- dependent approach is 
to construct a completely local theory, for which local re- 
flection and transmission amplitudes associated with any 
given trajectory, as it encounters a given step k, depend 
only on the properties of the fc'th step (i.e. on Xk, and on 
the p or V values to the immediate left and right oi Xk). 



10 




1500 



time (a.u.) 



FIG. 3: Bipolar trajectory plot for the square barrier prob- 
lem with E > Vo, as described in Sec. IIVBI One trajectory 
in five is indicated in the figure. The black/gray solid lines 
indicate positive/negative LM trajectories, respectively. The 
open circles represent recombination points. The dashed lines 
denote the two barrier edges, xi = and X2 = 1. 



In fact, from the point of view of the given trajectory, it 
must be immaterial whether the potential contains other 
steps or not — implying that the correct local relations 
for the spawned trajectories, if they exist at all, must be 
exactly those already specified in Eqs. ([T7]) and p^ . 

How is it possible that for stationary wavefunctions, 
whose time evolution is presumably trivial, an inherently 
global problem can be converted to a local one, simply by 
switching from a time-independent to a time-dependent 
perspective? This is because of the bipolar decompo- 
sition, which provides each step with not one, but two 
sets of incident trajectories, one from the left, and one 
from the right. When there are multiple steps, not only 
does this result in a non-trivial superposition for the re- 
sultant locally reflected and transmitted waves, but the 
trajectories themselves are subject to multiple spawnings, 
which effectively enable them to traverse back and forth 
over the same regions of configuration space an arbitrary 
number of times (Fig. [3]). This crucial feature ultimately 
gives rise to the rich global scattering behavior observed 
even in two-step systems. However, it is wholly missed 
by any time-independent treatment, even a bipolar one, 
which can only summarize the net superposition of all 
left-traveling and right-traveling waves. 

We now discuss how the local time-dependent theory 
described above gives rise to the correct stationary so- 
lutions, which is readily understood by invoking the ray 
optics description introduced earlier. For simplicity and 
definiteness, we consider only the square potential case, 
which is optically analogous to say, a single pane of glass 
surrounded by vacuum. If a single step gives rise to a sin- 
gle reflection, then two steps, like a pair of mirrors, results 
in an infinite number of reflections. The same is true of a 
pane of glass, within which a single beam of light will be 
reflected back and forth at the edges an arbitrary number 
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FIG. 4: Seven snapshots of the superposition wavefunction, 
'i/{x,t), for the E > Vq square barrier problem, as computed 
using the numerical algorithm of Sec. lIIII The shaded box rep- 
resents the barrier region. All units are atomic. The evolving 
discontinuities found at intermediate times in these curves de- 
note wavefronts for the newly created reflected/transmitted 
components. Over time, the magnitudes of these disconti- 
nuities (i.e. corrections) become arbitrarily small, signifying 
that numerical convergence to the left-incident stationary so- 
lution has been achieved. 



of times. Of course, these reflections are not perfect; a 
portion of the incident flux always escapes as transmis- 
sion into the surrounding vacuum. Consequently, each 
successive internal reflection is exponentially damped, in 
accord with Eq. (fTT)) . 

If the globally incident wave is incoming from the left, 
then at xi, there are two contributions to b+- One con- 
tribution is the portion of the left-incident '^a+ wave that 
is locally transmitted through the flrst step. Apart from 
a phase factor, the resultant value would be given 
by Eq. (jlSp if this were the only contribution. However, 
there is also a contribution that arises from the locally 
reflected part of the right-incident wave, \E'b-- This con- 
tribution is zero for a single step system, but of course 
non-zero in the multiple step case. Although the sec- 
ond contributing wave is right-incident, we can still use 
Eqs. (fTT)) and (fT5)) to compute the contribution to ^'b-i-, 
as discussed in Sec. Ill DT] For the second, k = I = 2 step 
at X2, there are only left-incident waves; consequently, 
^(7+ and '^B- are obtained from a single source each, 
i.e. (Fig.©. 

The above description refers to the stationary state re- 
sult, obtained by our gedankenexperiment in the large 
time limit only. In practice this result would be achieved 
in stages. As in the previous examples, we imagine that 
at time t = 0, one retains only those trajectories for 
which X < XL < xi. This one-sided trajectory restric- 
tion is somewhat analogous to continuous wave cavity 
ring-down spectroscopy.^^ When the wavefront flrst hits 
the flrst interface at a; = cci, there is partial reflection 
and transmission, exactly identical to what would hap- 
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pen for a single step system. The reflected wavefront 
propagates beyond the left edge of the region of interest 
at X = XL, and for some time, the reflected amplitude 
passing through this left edge is constant. The initially 
transmitted wavefront eventually reaches the second step 
a.t X — X2 (i.e. the far side of the pane of glass), lead- 
ing to a second transmission into the C region, and a 
second reflection back through the B region. Eventu- 
ally, the second transmitted wavefront reaches the right 
edge of interest at a; = xr, after which the transmitted 
amplitude remains constant for some time. 

Neither the globally transmitted nor reflected ampli- 
tudes for the times indicated above, as determined via 
monitors at x = Xl and x — xr respectively, are cor- 
rect. However, we have not yet described the steady 
state solution. To do so requires an accounting of the 
second reflected wavefront, which eventually reaches the 
X — Xl step again, this time incident from the right. 
The resultant locally transmitted wave becomes an in- 
stantaneous second contribution to ^'a-, and the locally 
reflected wave plays the same role for ^b+- These new 
contributions give rise to discontinuities in these waves, 
that subsequently propagate to the left and right, respec- 
tively (Fig[4]). The new ^ a- wave discontinuity eventu- 
ally reaches x — xl, where it is recorded by the monitor, 
giving rise to a sudden change in the reflection probabil- 
ity value. 

The vj/ discontinuity propagates to the second step, 
where it spawns new discontinuities in '^c+ and ^'s-. 
The former constitutes the border between first- and 
second-order transmitted waves, registered at sufficiently 
later time by the monitor at a; = xp^. The latter, second- 
order ^ B- wave heads back towards the first step, to give 
rise to third-order waves, with commensurate disconti- 
nuities, etc. In principle, this process continues indefi- 
nitely, resulting over time in global transmitted and re- 
flected waves of arbitrarily high order. However, Eq. p7|) 
and the relation Picfl + ^tians = 1 ensure that the result 
converges to a stationary solution exponentially quickly. 
Moreover, since C_ is necessarily zero throughout this 
process, it is clear that the stationary state that is con- 
verged to is indeed the one corresponding to the desired 
boundary condition of a globally incident wave that is 
incoming from the left. 

Note that in an actual optical system as described 
above, the spatial dimensionality is three rather than 
one, and the incident wave would usually be taken at 
some angle to the normal. If in addition, the beam has 
a finite width, then one would observe separate reflected 
beams for each order, of exponentially decreasing bright- 
ness. The one-dimensional quantum case, however, is 
analogous to a normal incident beam, for which all or- 
ders of reflection are superposed. In addition to provid- 
ing a pedagogical understanding of the dynamics that is 
very much analogous to the optical example provided, the 
picture above also suggests a practical numerical method 
that may be used to obtain stationary scattering states 
of any desired boundary condition (via superposition of 



globally left- and right-incident wave solutions, obtained 
independently) . 

Note that the "wavepacket dynamics" version of the 
ray optics analogy may also be applied. In this case, 
the resultant initial square wavepacket is somewhat rem- 
iniscent of pulsed wave cavity ring-down spectroscopyi^ 
Once the wavepacket has penetrated the middle region B 
(i.e. the pane of glass), it reflects back and forth between 
the two edges, with each reflection giving rise to a left- or 
right-propagating outgoing square wavepacket in region 
A or C, and a temporary interference pattern in region 
B. The amplitude of the central wavepacket dissipates 
exponentially in time. All of this complicated behavior 
is indeed qualitatively observed in actual wavepacket dy- 
namics for such systems, but in the present context, is 
reconstructed entirely from a single stationary state. 



III. NUMERICAL DETAILS 

In this section, we discuss several remaining issues 
pertaining to the numerical methods used to generate 
and propagate the various bipolar component waves, for 
the examples discussed in Sees. IIIDI and IIVI For all 
of these examples, the numerical algorithm used corre- 
sponds to the gedankenexperiment with one-sided trun- 
cation, i.e. to continuous wave cavity ring-down. In 
essence, this consists of just two basic operations: (1) 
each piece-wise bipolar component of the wavefunction 
[i.e. '^A±{x), '^B±{x), etc.] is independently propagated 
in time over its appropriate region of space, using the 
standard QHEMs and QTMs; (2) whenever a trajectory 
reaches a turning point, it is immediately deleted, and re- 
placed with two new trajectories, spawned in the appro- 
priate locally transmitted and reflected component LMs. 

The first operation above, i.e. QTM propagation of 
the wavefunction components, is very straightforward. 
Note that for simplicity, we have throughout this paper 
used time- independent expressions for '^{x) and its com- 
ponents, but in reality these evolve over time — even for 
stationary states, via s = ds{x,t)/dt = —E. We have 
therefore been rather lax in distinguishing Hamilton's 
principle function from Hamilton's characteristic func- 
tion, although from a trajectory standpoint, it is always 
the former that is implied. Since each component is sta- 
tionary in its own right, the time evolution of the hy- 
drodynamic fields is governed by the quantum station- 
ary Hamilton- Jacobi equation (QSHJE), rather than the 
QHJE. Moreover, the fact that the piecewise r is con- 
stant implies that the component quantum potentials are 
zero, resulting in classical HJE's and trajectories. These 
conclusions are trivially correct for the present paper, 
for which all components are plane waves; however, the 
arguments also extend to arbitrary continuous potential 
systems.**^ 

From a numerical perspective, the use of classical 
trajectories offers many advantages over a conventional 
QTM propagation. To begin with, the trajectories them- 
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selves are always smooth if V{x) is smooth, resulting in 
far fewer trajectories and larger time steps than would 
otherwise be the case. Even more importantly, however, 
since the quantum potential is not required, there is no 
need to compute on-the-fly numerical spatial derivatives 
of the local hydrodynamic fields. Consequently, for a 
given component, the trajectories are completely inde- 
pendent and need not communicate — again resulting in 
fewer of them. Indeed, it is possible to perform an essen- 
tially exact computation using only a single trajectory 
per wavefunction component. This feature is particularly 
important for the very frontmost trajectory of the initial 
ensemble, which for brief periods at later times, will (via 
spawning) come to be the only trajectory to occupy a 
given component LM. The subsequent evolution of these 
lone trajectories does not require the presence of nearby 
trajectories. 

The spawning of new trajectories, i.e. operation (2) 
above, also bears further discussion. In principle, this 
is always achieved via application of Eqs. (fT7|) and p8| . 
For the above barrier case, R and T are both real, ensur- 
ing the reality of r and s for the spawned trajectories — 
although in the case of a down step, i? < 0, resulting in 
a negative rj-ofl. This is in accord with the conventions 
discussed in paper I. However, in this paper, we find it 
numerically convenient to adopt the more usual r > 
convention. Thus, if Eq. ([TT]) yields a negative rj-cfl, it is 
replaced with — Ti-cfi, and tt is added to Srofl- A similar, 
but more complicated modification is also applied to the 
below-barrier trajectories, for which Eqs. pT]) and 
yield complex amplitudes. In this case, the phase shift is 
26, as discussed in Sec. Ill D"2] and the Appendix. 

For a single step system, the algorithm is now essen- 
tially complete. At the initial time t = 0, a variable 
number of particles (or synonymously, grid points) are 
distributed uniformly along the manifold, to the 

left of a; = xl- The extent of these points must be large 
enough that at the end of the propagation, there are still 
grid points that have not yet reached xl- The grid 
spacing is mostly arbitrary, but must be small enough 
that at sufficiently later times, there is always at least 
one trajectory per component LM. The propagation is 
considered complete when the reflected and transmitted 
wavefronts travel beyond xl and xr, respectively. 

For multiple step systems, the situation is similar, but 
somewhat more complex. The primary new feature is 
the recombination of wavefunction components arising 
from two sources, i.e. from two locally incident waves 
coming from opposite directions. The present algorithm 
would seem to yield two subcomponent wavefunctions 
for every component, each with its own set of trajecto- 
ries. If left "unchecked," this would lead to undesirable 
further multifurcations for higher orders/later times. A 
simple solution would be to propagate each subcompo- 
nent long enough that there is at least one trajectory 
for each, then extrapolate the corresponding subcom- 
ponent wavefunctions to a common position, where a 
new trajectory is constructed for the superposed com- 



ponent wavefunction, which is then propagated in lieu 
of the subcomponent trajectories. This requires dynam- 
ical fitting (see below), or at the very least, extrapo- 
lation. Although these numerical operations would be 
very stable in the present context, to rule these out al- 
together as sources of error in Sec. IIV[ we have adopted 
a much simpler approach — i.e. the grid spacing is cho- 
sen such that trajectories from the two component waves 
incident on a given step always arrive at the same time 
(Fig. [3]). The corresponding subcomponent wavefunc- 
tion values are then simply added together when form- 
ing the spawned trajectory. Adopting once again the 
r > convention for the superposed component wave, 
^1^-1-, the corresponding field values are then obtained via 
r = y¥|¥± and s = arctan [Im(\I'±)/Re(*±)]. 

As discussed in Sec. Ill multiple step systems allow 
for infinite reflections that perpetually modify |\E'(x,t)p, 
in principle for all time. In practice however, there is 
exponential convergence within the region of interest, 
XL ^ X < xr, so that one would not run the calcula- 
tion indeflnitely, but only until the desired accuracy is 
reached. Accurate "error bars" on the computed global 
Prcfl and Ptrans valucs are conveniently provided by the 
magnitudes of the most recent discontinuous jumps as 
recorded by the monitors at x^ and xr. Note that the 
number of digits of accuracy scales only linearly with 
propagation time. However, the rate of convergence de- 
pends on the energy value. Near the barrier height, in 
particular, convergence may take quite a long time, as 
the exponent is close to zero. For all other energies, only 
a few "cycles" should be required, depending on the level 
of accuracy desired. 



If in addition to reflection and transmission probabili- 
ties, the actual stationary solution over the region of in- 
terest is also desired, then it is necessary to reconstruct 
"^{x). This is obtained from the final grid, after the prop- 
agation is finished, using a multiple step generalization of 
Eq. (jl4p . The first step is to reconstruct the component 
wavefunctions '^a±{x), '^b±{x), etc., via interpolation 
or fitting of the hydrodynamic field values from the cor- 
responding dynamical grid points onto a much finer com- 
mon grid (used e.g. for plotting purposes). The second 
step is to linearly superpose the ± components onto the 
plotting grid, and to assemble the pieces together over 
the coordinate range of interest. For the discontinuous 
systems considered here, the number of dynamical grid 
points per component can be as small as one — i.e. much 
smaller, even, than the number of wavelengths! To our 
knowledge, such performance has never been achieved 
previously by a QTM; however, it does require that the 
plotting grid be much finer than the dynamical grid, e.g. 
at least several points per wavelength, in order to ade- 
quately represent the interference fringes of the the su- 
perposed solution, ^'(a;). 



13 



IV. RESULTS 

In this section, we apply the numerical algorithm pre- 
viously described to three different applications: the up- 
step potential, the square barrier, and the square well. 



A. Up-step Potential 

The first system considered is the up-step potential, 
i.e. Eq. pT|) with Vq > 0. Since there are no multiple 
reflections, this is in principle a trivial application for 
the current algorithm; it therefore serves as a useful nu- 
merical test. Both above barrier (Sec. HTdTI) and below 
barrier (Sec. Ill D 2p energies are considered. We choose 
molecular-like values for the constants, i.e. Vq = 0.009 
hartree, and m — 2000 a.u. The left and right edges of 
the region of interest are taken to he xl = — 1.0 a.u. and 
xr = 1.0 a.u., respectively. At the initial time, t — 0, 51 
trajectory grid points are distributed uniformly over the 
interval —4 < x < — 1 (grid spacing of 0.06 a.u.). This 
number is far greater than what would be needed for dy- 
namical purposes, but is chosen so as to avoid construc- 
tion of a separate plotting grid (Sec. IIIip . The hydrody- 
namic field functions for the initial "^A+ix) wavepacket 
over the above interval are taken to be r{x) — 1 a.u.^^ 
and s{x) — \/2mEx. 

For the above barrier calculation, the energy E — 
2Vq = 0.018 hartree was used. The trajectory propa- 
gation and termination were performed exactly as de- 
scribed in Sec. IIIII The real and imaginary parts of all 
three resultant wavefunction components [i.e. 'i'A±{x) 
and "^B+i^)] at the final time, t = 550 a.u. are pre- 
sented in Fig. [TJ All three components exhibit the de- 
sired plane wave behavior, e.g. no interference is evident 
within a given component. The resultant "^{x) does ex- 
hibit interference in the A region, however, arising from 
the superposition of '^a+ and '^'a-- 

For the below barrier calculation, the system was given 
an energy equal to one half of the barrier height, i.e. 
S = Vb/2 0.0045. As per the discussion in Sees. |IIDJ 
and lllli tunneling into the forbidden region B is achieved, 
not through a quantum potential, but via analytic con- 
tinuation. At sufficiently large time {t — 1100 a.u.) the 
final wavefunction is reconstructed from the components, 
i.e. ^^(x) = *A+(a;) + ^^-(a;), and ^'^(a;) = 1'B+(a;). 
The real and imaginary parts of the reconstructed wave- 
function are presented in Fig. [21 In the figure, squares 
and circles denote the numerical results obtained via the 
present algorithm, whereas the solid and dashed lines rep- 
resent the well-known analytic solutions. The agreement 
is essentially exact. Note that the real and imaginary 
parts are in phase throughout the coordinate range — i.e., 
S{x) is a constant, so apart from a phase factor, ^(x) is 
real. Note also that the tunneling region exhibits the 
desired exponential decay. 



B. Square Barrier 

The second system considered is the square barrier. 
This is a two-step potential {I — 2), with Va = Vc = 
0, and Vb = Vb > 0. The two steps comprise the left 
and right edges of the barrier, at xi = and X2 = w, 
respectively. The constants are chosen as follows: Vb = 
0.018 hartree; m = 2000 a.u.; w — 1] xl — —1 a.u., xr — 
2 a.u. Initially, 75 trajectory grid points are distributed 
uniformly over the interval — 5 < a; < — 1 (grid spacing of 
0.05 a.u.), which again, is far more than are dynamically 
required. The same initial hydrodynamic field functions 
are used as in Sec. IIV Al 

Both above barrier {E > Vq) and below barrier {E < 
Vq) energies are considered. For the above barrier case, 
E = 2Vo — 0.036 hartree. Once again, the trajectory 
propagation and termination were performed exactly as 
described in Sec. IIIII In order to converge Ptrans to 10~*, 
a propagation time of 3000 a.u. was required. This corre- 
sponds to 3 complete cycles, i.e. a 3rd-order calculation. 
Figure [3] is a plot of the quantum trajectories for this 
calculation, in which every fifth trajectory for each of 
the five component wavefunctions is indicated. Trajec- 
tory spawning at the two steps is very clearly evident, 
as is recombination of pairs of incident waves (indicated 
by circles). On the whole, this figure demonstrates all 
of the anticipated analogues with ray optics, i.e. parallel 
trajectories, reflection and refraction. 

In Fig. 31 the time evolution of the superposition state 
^'(a;, t) is represented, via snapshots of the real and imag- 
inary parts at seven different times. At t = a.u., the 
incident wavefront is located at x = —1 a.u. By t = 250 
a.u., the wavefront has spawned '^a-{x) and ^b+{x) 
trajectories; the former gives rise to the kink (really a 
discontinuity) somewhat to the left of the first step. By 
t = 550 a.u. and t = 650 a.u., the ^'^-(a;) wavefront has 
moved outside the region of interest, though the 5*5+ (.x) 
wavefront has not quite reached the second step. After 
it does so, two new wavefronts are propagated along the 
^'c+(x) and ^'B_(a:) LMs (e.g. t = 800 a.u.), the for- 
mer of which propagates beyond the region of interest by 
t = 900 a.u. Subsequent discontinuity magnitudes be- 
come exponentially smaller, so that at sufficiently large 
time (i.e. t = 2000 a.u.), the resultant ^'(a;) has con- 
verged to the correct stationary solution. 

For the below barrier case, E = Vo/2 = 0.009 hartree, 
XL — —0.5 a.u., and the other parameters are as above 
except w = 0.5 a.u. Within the barrier, there are in prin- 
ciple an arbitrary number of reflections back and forth as 
before. However, substantial amplitude loss occurs due 
to tunneling, in addition to partial reflection, as a result 
of which fewer cycles are required in order to achieve 
the same 10~^ level of convergence {t = 1400 a.u., or 
2 cycles). Fig. |S1 indicates how the tunneling dynamics 
is achieved. After the wavefront hits the flrst step, the 
transmitted "^3+ wave is propagated along the imaginary 
axis (iy), until the point y = w is reached. When this 
occurs, it is necessary to analytically continue down 
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FIG. 5: Schematic of the algorithm used to propagate tra- 
jectories into the classically forbidden region of the E < Vo 
square barrier problem, as discussed in Sees. IIIlland llVB] 



to the real axis, in order to compute amplitudes for the 
new and ^b-{x) trajectories that are spawned 

at the second step. The latter component propagates in 
the (negative) imaginary direction, along w — iy' , until 
y' = w, at which point analytic continuation is once more 
applied (this time for the first step), and the pattern re- 
peated. 

Seven snapshots of the superposition density, 
|4'(x,t)p, are displayed in Fig. [6] The initial density — 
equal to just the |\['^+(a;)p density — is uniform. After 
the wavefront encounters the first step, a reflected 
"ifA-ix) emerges, giving rise to clearly evident inter- 
ference in the region A. The ^b{x) = ^B+ix) wave 
is at this stage perfectly exponentially damped. Upon 
encountering the second step, a second contribution, 
^B-ix) emerges; however, this first-order correction is 
already extremely small, owing to the large amount of 
tunneling that has occurred. The global transmitted 
wave, \E'c-i-(a:), though small, is clearly seen to have 
uniform density. 

In addition to the two detailed trajectory calculations 
described above, we computed Prcfl and Ptrans for a large 
range of w and E values, so as to fully explore (without 
loss of generality) the entire range of the square barrier 
problem. The numerical results are presented, and com- 
pared with known analytical valueSf^ in Fig. [71 Two 
aspects of this study bear comment. First, for all w and 
E values considered, the computed Profl and Ptrans values 
agree with the exact values to within an error compara- 
ble to that predicted by the level of numerical conver- 
gence. In particular, the oscillatory energy dependence 
is perfectly reproduced. Second, the closer the barrier 
peak is approached from either above or below in energy, 
the longer the time required to achieve a given level of 
convergence, as predicted in Sec. IIIIl In particular, for 
the calculations closest to the barrier peak, 5-7 cycles 
were required in order to approximately maintain a 10"'' 
convergence of the transmission probability. Although a 
greater number of particles are required in this case, this 
poses no great limitation in practice, since one would pre- 
sumably never require a calculation precisely at the peak 





FIG. 6: Seven snapshots of the superposition density, 
|*(a:,i)P, for the E < Vo square barrier problem, as com- 
puted using the numerical algorithm of Sec. IIIIl The shaded 
box represents the barrier region. All units are atomic. Note 
that interference manifests in the incident (left) region only af- 
ter some incident trajectories have struck the left barrier edge, 
causing reflected trajectories to be created (i.e. just prior to 
t — 195, initially). The most advanced reflected trajectory de- 
fines the reflected wavefront, manifesting as the left-moving 
discontinuity, e.g. at t = 195 and t — 287. 
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FIG. 7: Transmission and reflection probabilities as a func- 
tion of energy, for square barrier potentials of three different 
widths, w, as discussed in Sec. IIVBI Solid lines denote an- 
alytical results; open/closed circles denote numerical results, 
as obtained via algorithm of Sec. IIIIl The vertical dot-dashed 
lines represent the barrier height, i.e. E = Vq. 



energy. 
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FIG. 8; Transmission and reflection probabilities as a func- 
tion of energy, for square well potentials of three different 
widths, w, as discussed in Sec. IIV CI Sohd hues denote an- 
alytical results; open/closed circles denote numerical results, 
as obtained via algorithm of Sec. IIIII 



C. Square Well 

As the final system, we consider the square-well po- 
tential, i.e. the square barrier but with Vq < 0. In the 
scattering state context, there is no tunneling for this 
system, but in other respects it resembles the square 
barrier. From an optics point of view, the square well 
corresponds to a central medium with larger index of re- 
fraction than its surroundings, whereas the square barrier 
corresponds to a smaller index of refraction, giving rise 
to the possibility of total internal reflection (i.e. tunnel- 
ing). The parameters are as in Sec. IIVBI except that 
Vq = 0.009 hartree, and three different w values are con- 
sidered: w = 2 a.u., w = A a.u., and it; = 16 a.u. 

As the time evolution and trajectory pictures are sim- 
ilar to those of the previous sections, we focus only on 
the Profl and Ptrans calculations, for which once again, a 
large range of energies was considered (0.0005 < i? < 0.2 
hartree). The number of initial trajectories ranged from 
50 to 200, for the highest to the lowest energies, respec- 
tively. The computed transmission/reflection probabili- 
ties were again converged to 10^*. The energy-resolved 
reflection and transmission probabilities are presented in 
Fig.m 

As in the square barrier case, excellent agreement is 
achieved with the exact analytical results, i.e. on the or- 
der of the level of convergence. This is true despite the 
fact that the square well energy curves are decidedly more 
oscillatory than the square barrier curves — particularly 
for wide barriers, for which the E dependence is very 
sensitive indeed. One particularly important feature ex- 
hibited by the exact curves is the so-called Ramsauer- 



Townsend effect;^ i.e. the phenomenon of 100% trans- 
mission and zero reflection, even at very low energies. 
This occurs when sin(2fcBw) — 0, and may be regarded 
as a purely quantum mechanical resonance phenomenon. 
Yet it is reproduced perfectly here in the bipolar decom- 
position, using classical trajectories. Indeed, the bipolar 
picture provides an interesting physical explanation, i.e. 
the right-incident and left-incident waves of the first step 
give rise to spawned contributions to "^A-ix) that ex- 
actly cancel each other out via destructive interference. 



V. SUMMARY AND CONCLUSIONS 

As described in paper I, the Schrodinger equation is 
linear, yet the equivalent QHEM — obtained via substitu- 
tion of the Madelung-Bohm ansatz into the Schrodinger 
equation — are not. This aspect of Bohmian mechanics 
suggests that it can be beneficial, both from a pedagog- 
ical and a computational perspective, to apply a suit- 
able bifurcation (or "multifurcation" ) to the wavefunc- 
tion prior to applying the QHEM. Indeed, following the 
paper I CPWM bipolar decomposition for ID stationary 
bound statesj^^ the quantum trajectories become more 
well-behaved and classical-like, in precisely the limit in 
which there are more nodes, and the usual unipolar calcu- 
lation breaks down. Moreover, the resultant component 
LMs admit a natural physical interpretation in terms of 
the corresponding semiclassical LMs. 

In the generalization to the stationary scattering states 
considered here, a somewhat different bipolar decompo- 
sition is found to be required. The new decomposition is 
still unique, at least for discontinuous potentials. How- 
ever, the resultant components ^±{x) are no longer so- 
lutions to the Schrodinger equation in their own right, 
as a result of which their time evolution is coupled. The 
new scheme — though fundamentally different from the 
old one — nevertheless bears a correspondence to a mod- 
ified version of semiclassical theory appropriate for scat- 
tering systems. Curiously, this semiclassical modification 
is not simply a higher order treatment in if it were, 
the corresponding exact quantum modification consid- 
ered here would not exist. 

In any event, the new decomposition also gives rise 
to its own physical interpretation, specific to the scat- 
tering context. In particular, for right-incident bound- 
ary conditions, the left and right asymptotes of ^'-i-(a;) 
respectively represent incident and transmitted waves, 
whereas the left asymptote of ^'_(a;) represents the re- 
fiected wave [the right asymptote of '^-{x) approaches 
zero]. That these conceptually useful asymptotic bipo- 
lar assignments — found even in the most elementary 
treatments of scattering — may be extended throughout 
configuration space [even for continuous potentials (pa- 
per HI)] represents an important leap forward, especially 
for QTMs. 

Another pedagogically and numerically useful devel- 
opment from this approach is the inherently time- 
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dependent ray optics interpretation that naturally arises, 
particularly in the discontinuous potential context. The 
ray optics approach is anticipated to be a relevant guid- 
ing force in subsequent generalizations of the present 
methodology, i.e. to continuous, multidimensional po- 
tential systems, and — as per the discussion in Sees. Ill 
and HID 11 — for non-stationary wavepacket dynamics. 
This approach also provides a much simpler trajectory- 
based explanation of scattering terminology as applied 
in a stationary context — e.g. why the "reflected wave" is 
so-called, despite being present from the earliest times — 
than those traditionally usedi^ 

The discontinuous potential applications considered 
in this paper — hard wall (Sec. Ill C 2[) . step potential 
(Sees. Ill D II HID 21 and IIV A|) . and square barrier/weU 
(Sees. IiTd31 IIV B| and II V C]) — are significant for several 
reasons. To begin with, these are the first discontinuous 
applications of a genuine QTM calculation that have ever 
been performed, to the authors' knowledge. The singular 
derivatives associated with discontinuous potential func- 
tions would wreak havoc with standard numerical differ- 
entiation routines. Second, the use of a time- dependent 
method for stationary, or time-independent, applications, 
is also significant. Ordinarily, the time dependence of 
stationary states is regarded as trivial. In the present 
context, this is true in a sense for the hard wall and step 
potential systems, because the correct answer is "built 
in" the method itself. For multiple step systems, how- 
ever, the dynamical truncated wave approach, i.e. the 
gedankenexperimcnt introduced in Sec. Ill C"2l and further 
developed in later sections, yields decidedly nontrivial re- 
sults. 



In particular, the algorithm uses only single step scat- 
tering coefficients to obtain global scattering quantities 
for multiple step systems. In effect, the time depen- 
dent nature of this approach allows computation of global 
properties using a completely local method. Not only 
were exact quantum results obtained for a full range of 
system parameters, but the numerical resources neces- 
sary to achieve this — i.e. the number of trajectories and 
time steps — were decidedly minimal. Indeed, the algo- 
rithm lives up to the promise made in paper I, of per- 
forming an accurate quantum calculation with fewer tra- 
jectories than nodes — a prospect virtually unheard of in 
a unipolar context. 

In future publications, we will naturally attempt to 
generalize the methodology described here and in pa- 
per III, for the type of multidimensional time-dependent 
wavepacket dynamics relevant to chemical physics appli- 
cations. In this context, the scattering version of the 
CPWM decomposition developed here is an absolutely 
essential first step, as reactive scattering is the underpin- 
ning of all chemical reactions. Additional discussion and 
motivation will be provided in paper III. 



ACKNOWLEDGEMENTS 

This work was supported by awards from The Welch 
Foundation (D-1523) and Research Corporation. The au- 
thors would like to acknowledge Robert E. Wyatt and 
Eric R. Bittner for many stimulating discussions. David 
J. Tannor and John C. TuUy are also acknowledged. 



Appendix: Bipolar decomposition of semi-bound 
states 



As discussed in Sees. HlBl III C 21 and III D 21 semi- 
bound stationary states in ID are bounded on one side 
only, as a result of which they are real- valued and singly- 
degenerate, like bound states. Consequently, they are 
amenable to the CPWM bipolar decomposition scheme 
introduced in paper I. In this appendix, we apply this de- 
composition to two semi-bound systems: the hard wall 
system, and the below-barrier up-step system. 



A. Hard wall system 

From Ref. [l^, the most general bipolar decomposi- 
tion of a hard wall stationary state — corresponding to 
Sec. Ill A"2l condition (1) only — is found to satisfy 



cot [s{x)/h] 



/mF 



— cot{kx) 



B 



(19) 



where s(x) 



-s_(a;), and ^+(2;) = r_(a;) is 



obtained from s{x) via Eq. ([7|) (without "sc" subscripts). 
The arbitrary parameters F and B are the invariant 
flux and median action parameters associated with con- 
ditions (2) and (3), respectively,^^ although the definition 
of F has been changed slightly to account for the scat- 
tering normalization convention, r^{x) — 1. Note that 
only the scmiclassical values for these parameters yields 
a solution that satisfies the correspondence principle in 
the large action (i.e. k) limit. In particular, the choice 
B — and F = Tik/m yields the desired scmiclassical 
result, s{x) = hkx; all other choices exhibit undesirable 
oscillatory behavior in r±{x), s±{x), and q±{x). 



B. Up-step system 

For the hard wall system considered above — which is 
just the special case of the up-step potential in the limit 
Vq — > 00 — exact agreement is achieved between semiclas- 
sical and quantum LM's in the x < region. This is the 
only region of interest for the hard wall system; however, 
for finite Vq values — i.e. for general below-barrier up-step 
stationary states — there is of course also tunneling into 
the forbidden region, which must be accounted for. The 
paper I bipolar decomposition therefore results in LM's 
that span the entire coordinate range — 00 < a; < 00. 
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These LMs are given by the following analytical expres- 
sion: 



p{x) 



hk 



2hK.e^'^'' sin 25 



l+e-*»^-2e^"=" cos 2(5 

where 5 is given by 



for X < 0; 
for a; > 0. 




and 



K = \/2m(Vb - E). 



(20) 



(21) 



(22) 



The p{x) LM function of Eq. (|20p is continuous every- 
where, including at the potential discontinuity at x = 0. 
In the A region, it agrees exactly with the semiclassical 
solution; in the B region, it decays exponentially to zero. 

The paper I approach thus yields a very natural way to 
extend trajectories into the tunneling region. Note that 



the quantum potential in this region is not zero; indeed, 
it exhibits a discontinuity at a; = that exactly balances 
that of V{x) itself, so that the bipolar modified potential 
is continuous across the step. Unlike the above-barrier 
case, the paper I solution docs not manifest oscillatory 
behavior in the large action limit, and so this approach 
would at first glance appear to be ideal. There are two 
reasons, however, why it is not pursued here. The first 
reason is that r(x) diverges asymptotically as a; oo, 
which according to preliminary numerical investigations, 
appears to lead to numerical instabilities for completely 
QTM-based propagation schemes. Second, if the barrier 
were to fall off again at larger x values, so that the tun- 
neling region were finite, then the asymptotic behavior 
would be once again undesirably oscillatory. This would 
be the case, for example, for the below-barrier energies 
of the square barrier system of Sec. IIVBI 
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